%****************************************************************
%** This program is specifically written to replicate Xtreg in **
%** Hansford and Gomez's (2010) APSR article.                  **
%** Created by Min-Hua Huang, Jan 25, 2012                     **
%** Model: Maximum Likelihhod by Fixing at Minimum, with       **
%**        Constrained Optimization, Conducting Correction on  **
%**        Location Estimates.                                 **
%****************************************************************
clear; clc; format long;
load HansfordGomez_Data.mat;
load recordn.mat;
n=size(HansfordGomez_Data(:,27),1);
code=HansfordGomez_Data(:,2);
obs=(1:n)';
y1=[HansfordGomez_Data(:,27) code obs];
y1=sortrows(y1,2); m=size(recordn,1)-1; 
n=size(y1,1); yir0=ones(m,2); yimean=ones(m,2);
order=y1(:,2:3);
for i=1:1:m
    first=recordn(i); 
    last=recordn(i+1)-1;
    y=y1(first:last,1);
    r0=[mean(y);std(y)]; a=0; b=100; runmax=101; tolv=10^-4;
    fs1=trmcod(y,r0,a,b,runmax,tolv);    
    yir0(i,:)=fs1';
end
r1=[mean(yir0(:,1));std(yir0(:,1))];
fs2=trmcod(yir0(:,1),r1,a,b,runmax,tolv);
y2r0=fs2';
newy1=ones(n,1); 
for i=1:1:m
    first=recordn(i); 
    last=recordn(i+1)-1;
    newy1(first:last,1)=yir0(i,1);
end
temp1=[newy1 order];
temp2=sortrows(temp1,3);
y2=temp2(:,1);
grandmeany=y2r0(1);
y1=sortrows(y1,3);
% *********************************************************************
x0=[HansfordGomez_Data(:,26) HansfordGomez_Data(:,3) HansfordGomez_Data(:,25) HansfordGomez_Data(:,29) HansfordGomez_Data(:,11:23)];
grandmeanx0=mean(x0);
obs=(1:n)';
replicate=[HansfordGomez_Data(:,27) x0 code obs];
replicate=sortrows(replicate,19);
dmean=dmean1(replicate,19);
mean4a=[(replicate(:,1:18)-dmean) replicate(:,19:20)];
mean4a=sortrows(mean4a,20);
wfmin=ones(n,2); wf1min=ones(n,2); 
for i=1:1:2
    wfmin(:,i)=mean4a(:,i+1)-min(mean4a(:,i+1)); 
end
cr1=HansfordGomez_Data(:,10).*wfmin(:,2);
cr2=wfmin(:,1).*wfmin(:,2);
afixx=[ones(n,1) wfmin cr1 cr2 HansfordGomez_Data(:,11:23)];
y=y1(:,1)-y2;
ols_afix=(afixx'*afixx)\(afixx'*y);
%**** Set up initial values for constrained optimization ****
sigma=sqrt((y-afixx*ols_afix)'*(y-afixx*ols_afix)/(n-m-1));
x=afixx; 
r0=[ols_afix;sigma];
runmax=31;
dmy=13;
a=0;
b=100;
tolv=10^-4;
y1=y;
frac=1;
tau=50;
%**** Executing maximim likelihood constrained optimization ****
%*****************  Without bias correcttion *******************
fs=ML2(y1,x,r0,a,b,runmax,dmy,tolv,tau,y2r0,yir0,frac);
m1=size(x,2)-1;
ml2b=fs(1:m1+1)';
ml2sigma=fs(m1+2);
ml2se=robuster(y1,afixx,[ml2b;ml2sigma],code);
ml2t=ml2b./ml2se;
hr=boundary1(y2r0,yir0,frac);
MLCO2ymax=max(x*ml2b+y2);
MLCO2ymin=min(x*ml2b+y2);
xmax=zeros(m1+1,1); xmin=zeros(m1+1,1);
   for i=1:1:m1+1
       xmax(i)=max(x(:,i));
       xmin(i)=min(x(:,i));
   end    
xlimit=[xmax xmin];
v=vsign(ml2b);
MLCO2ci=cpar(fs(1:m1+2)',a-hr(2),b-hr(1),x,xlimit,v,dmy);
MLCO2bv=sum(MLCO2ci>tolv);
MLCO2predict=range1(ml2b,hr(2),hr(1),x,xlimit,v,dmy,y2);
fprintf('\nReplication of the column 1 in Table 1, Hansford and Gomez (2010)\n')
fprintf('\nMaximum Likelihod Estimation by Fixing at Minimum with Constrained Optimization\n')
fprintf('This Model Performs Bias Correction in the Demeaning Operation\n')
fprintf('                                 beta      se       t             \n');
fprintf('Partisan Composition           % 7.3f  % 7.3f  % 7.3f\n',ml2b(2),ml2se(2),ml2t(2));
fprintf('Turnout                        % 7.3f  % 7.3f  % 7.3f\n',ml2b(3),ml2se(3),ml2t(3));
fprintf('Turnout*GOP Incumbent          % 7.3f  % 7.3f  % 7.3f\n',ml2b(4),ml2se(4),ml2t(4));
fprintf('Turnout*Partisan Composition   % 7.3f  % 7.3f  % 7.3f\n',ml2b(5),ml2se(5),ml2t(5));
fprintf('Yr52                           % 7.3f  % 7.3f  % 7.3f\n',ml2b(6),ml2se(6),ml2t(6));
fprintf('Yr56                           % 7.3f  % 7.3f  % 7.3f\n',ml2b(7),ml2se(7),ml2t(7));
fprintf('Yr60                           % 7.3f  % 7.3f  % 7.3f\n',ml2b(8),ml2se(8),ml2t(8));
fprintf('Yr64                           % 7.3f  % 7.3f  % 7.3f\n',ml2b(9),ml2se(9),ml2t(9));
fprintf('Yr68                           % 7.3f  % 7.3f  % 7.3f\n',ml2b(10),ml2se(10),ml2t(10));
fprintf('Yr72                           % 7.3f  % 7.3f  % 7.3f\n',ml2b(11),ml2se(11),ml2t(11));
fprintf('Yr76                           % 7.3f  % 7.3f  % 7.3f\n',ml2b(12),ml2se(12),ml2t(12));
fprintf('Yr80                           % 7.3f  % 7.3f  % 7.3f\n',ml2b(13),ml2se(13),ml2t(13));
fprintf('Yr84                           % 7.3f  % 7.3f  % 7.3f\n',ml2b(14),ml2se(14),ml2t(14));
fprintf('Yr88                           % 7.3f  % 7.3f  % 7.3f\n',ml2b(15),ml2se(15),ml2t(15));
fprintf('Yr92                           % 7.3f  % 7.3f  % 7.3f\n',ml2b(16),ml2se(16),ml2t(16));
fprintf('Yr96                           % 7.3f  % 7.3f  % 7.3f\n',ml2b(17),ml2se(17),ml2t(17));
fprintf('Yr2000                         % 7.3f  % 7.3f  % 7.3f\n',ml2b(18),ml2se(18),ml2t(18));
fprintf('Constant                       % 7.3f  % 7.3f  % 7.3f\n',ml2b(1),ml2se(1),ml2t(1));
fprintf('Sigma                          % 7.3f\n',ml2sigma);
fprintf('Max yhat                       % 7.3f\n',MLCO2ymax);
fprintf('Min yhat                       % 7.3f\n',MLCO2ymin);
fprintf('Theoretical Boundary Violations% 7.0f\n', MLCO2bv);
fprintf('N                              % 7.0f\n', n);
fprintf('\nPredicted Value by Varying        Min      Max\n');
fprintf('All Variables                  % 7.1f  % 7.1f\n',MLCO2predict(1,1),MLCO2predict(1,2));
fprintf('Constant                       % 7.1f  % 7.1f\n',MLCO2predict(2,1),MLCO2predict(2,2));
fprintf('Partisan Composition           % 7.1f  % 7.1f\n',MLCO2predict(3,1),MLCO2predict(3,2));
fprintf('Turnout                        % 7.1f  % 7.1f\n',MLCO2predict(4,1),MLCO2predict(4,2));
fprintf('Turnout*GOP Incumbent          % 7.1f  % 7.1f\n',MLCO2predict(5,1),MLCO2predict(5,2));
fprintf('Turnout*Partisan Composition   % 7.1f  % 7.1f\n',MLCO2predict(6,1),MLCO2predict(6,2));
fprintf('Yr52                           % 7.1f  % 7.1f\n',MLCO2predict(7,1),MLCO2predict(7,2));
fprintf('Yr56                           % 7.1f  % 7.1f\n',MLCO2predict(8,1),MLCO2predict(8,2));
fprintf('Yr60                           % 7.1f  % 7.1f\n',MLCO2predict(9,1),MLCO2predict(9,2));
fprintf('Yr64                           % 7.1f  % 7.1f\n',MLCO2predict(10,1),MLCO2predict(10,2));
fprintf('Yr68                           % 7.1f  % 7.1f\n',MLCO2predict(11,1),MLCO2predict(11,2));
fprintf('Yr72                           % 7.1f  % 7.1f\n',MLCO2predict(12,1),MLCO2predict(12,2));
fprintf('Yr76                           % 7.1f  % 7.1f\n',MLCO2predict(13,1),MLCO2predict(13,2));
fprintf('Yr80                           % 7.1f  % 7.1f\n',MLCO2predict(14,1),MLCO2predict(14,2));
fprintf('Yr84                           % 7.1f  % 7.1f\n',MLCO2predict(15,1),MLCO2predict(15,2));
fprintf('Yr88                           % 7.1f  % 7.1f\n',MLCO2predict(16,1),MLCO2predict(16,2));
fprintf('Yr92                           % 7.1f  % 7.1f\n',MLCO2predict(17,1),MLCO2predict(17,2));
fprintf('Yr96                           % 7.1f  % 7.1f\n',MLCO2predict(18,1),MLCO2predict(18,2));
fprintf('Yr2000                         % 7.1f  % 7.1f\n\n',MLCO2predict(19,1),MLCO2predict(19,2));
